# install.packages("pacman")
# install.packages("devtools")
# install.packages("magick")
# devtools::install_github("clessn/clessnize", force=T)
# devtools::install_github("clessn/hublotr", force=T)
# devtools::install_github("ellipse-science/sondr", force=T)
# devtools::install_github("clessn/clessnverse", force=T)
pacman::p_load(modelsummary, psych, tidyverse, clessnverse, haven, pandoc, knitr, kableExtra, marginaleffects, extrafont, patchwork, extrafont)

JaponW2 <- read.csv("_SharedFolder_article_technophobie/data/JapanW2cleaned.csv")
# 2023-07-04 14:59:13 transform from character to Date
JaponW2$StartDate <- as.Date(JaponW2$StartDate, format = "%Y-%m-%d %H:%M:%S")
summary(JaponW2$StartDate)

# 1. Explore and clean variables ####
## 1.1 AI dread
# table(JaponW2$q9_1, useNA = "always") # In the future, society will be dominated by artificial intelligence.
# table(JaponW2$q9_2, useNA = "always") # Artificial intelligence will surpass one day all aspects of human intelligence.
table(JaponW2$q9_3, useNA = "always") # I fear that artificial intelligence will attack and harm humans.
# table(JaponW2$q9_4, useNA = "always") # The use of artificial intelligence technologies in the workplace will lead to many job losses.
table(JaponW2$q9_5, useNA = "always") # I fear that artificial intelligence will ultimately replace all humans.
table(JaponW2$q9_6, useNA = "always") # Artificial intelligence poses a grave threat to our civilization.
table(JaponW2$q9_7, useNA = "always") # I fear that artificial intelligence will bring catastrophe such as the end of mankind.
## 1.2 AI controllability
# table(JaponW2$q10_1, useNA = "always") # I am concerned about artificial intelligence becoming too widespread in society.
table(JaponW2$q10_2, useNA = "always") # I believe that experts can control the impact of artificial intelligence.
table(JaponW2$q10_3, useNA = "always") # It is possible to accurately predict how artificial intelligence will affect different types of jobs.
table(JaponW2$q10_4, useNA = "always") # We can understand how much artificial intelligence will change society.
# table(JaponW2$q10_5, useNA = "always") # I fear that one day humans will be under the control of artificial intelligence.
table(JaponW2$q10_6, useNA = "always") # It is possible to limit the risks of artificial intelligence technology.
# table(JaponW2$q10_7, useNA = "always") # I am concerned that businesses will use artificial intelligence in ways they cannot control.
## 1.3 AI vignette (but unsure if dread or controllability)
table(JaponW2$q13_1, useNA = "always") # I am concerned A.I. will have negative effects on society and individuals.
table(JaponW2$q13_2, useNA = "always") # I think artificial intelligence technology will never fully match humans in performing tasks and making decisions.
table(JaponW2$q13_3, useNA = "always") # I think the government should oversee the development of new artificial intelligence technologies.
table(JaponW2$q13_4, useNA = "always") # Citizens should have a say in deciding how artificial intelligence is used in society.
## 1.4 Left-right ideology
table(JaponW2$q22, useNA = "always") # Left-Right Ideology
JaponW2$leftRight_selfReport_num <- (JaponW2$q22 - 1) / 10
JaponW2$leftRight_selfReport_num[JaponW2$leftRight_selfReport_num == 1.1] <- NA
table(JaponW2$leftRight_selfReport_num, useNA = "always")
JaponW2$leftRight_selfReport_num_no5 <- JaponW2$leftRight_selfReport_num
JaponW2$leftRight_selfReport_num_no5[JaponW2$leftRight_selfReport_num_no5 == 0.5] <- NA
table(JaponW2$leftRight_selfReport_num_no5, useNA = "always")
# table(JaponW2$q11_5, useNA = "always") # Japan should only admit skilled and educated workers as permanent immigrants.
table(JaponW2$q15_1, useNA = "always") # The government should increase taxes on companies’ profits.
table(JaponW2$q15_2, useNA = "always") # The government should increase benefits paid to people who have lost their jobs.
# table(JaponW2$q15_3, useNA = "always") # The government should increase restrictions on the number of temporary foreign workers allowed to work in Japan.
table(JaponW2$q15_4, useNA = "always") # The government should pay every citizen enough to cover the cost of basic necessities, regardless of whether they work or not.

## 1.5 Control variables ####
### 1.5.1 Gender ####
table(JaponW2$q1, useNA = "always")
JaponW2$ses_female_bin <- NA
JaponW2$ses_female_bin[JaponW2$q1 == 2] <- 1
JaponW2$ses_female_bin[JaponW2$q1 == 1] <- 0
table(JaponW2$ses_female_bin, useNA = "always")

### 1.5.2 Year of birth ####
table(JaponW2$q2, useNA = "always")
JaponW2$ses_age_num <- 2023 - JaponW2$q2
JaponW2$ses_age35l_bin <- NA
JaponW2$ses_age35l_bin[JaponW2$ses_age_num < 35] <- 1
JaponW2$ses_age35l_bin[JaponW2$ses_age_num >= 35] <- 0
JaponW2$ses_age3555_bin <- NA
JaponW2$ses_age3555_bin[JaponW2$ses_age_num >= 35 & JaponW2$ses_age_num < 55] <- 1
JaponW2$ses_age3555_bin[JaponW2$ses_age_num < 35 | JaponW2$ses_age_num >= 55] <- 0
JaponW2$ses_age55m_bin <- NA
JaponW2$ses_age55m_bin[JaponW2$ses_age_num >= 55] <- 1
JaponW2$ses_age55m_bin[JaponW2$ses_age_num < 55] <- 0
table(JaponW2$ses_age35l_bin, useNA = "always")
table(JaponW2$ses_age3555_bin, useNA = "always")
table(JaponW2$ses_age55m_bin, useNA = "always")

### 1.5.3 Immigrant status ####
# table(JaponW2$q4, useNA = "always")
# JaponW2$ses_personOrParentBornInJapan_bin <- NA
# JaponW2$ses_personOrParentBornInJapan_bin[JaponW2$q4 == 1] <- 1 # Person and parent born in Japan
# JaponW2$ses_personOrParentBornInJapan_bin[JaponW2$q4 == 2] <- 0 # Person or parent not born in Japan
# table(JaponW2$ses_personOrParentBornInJapan_bin, useNA = "always")
# UNUSABLE, TOO FEW CASES

### 1.5.4 Education ####
table(JaponW2$q5, useNA = "always")
JaponW2$ses_education_num <- JaponW2$q5
JaponW2$ses_education_num[JaponW2$ses_education_num == 7] <- NA
JaponW2$ses_educHighSchool_bin <- NA
JaponW2$ses_educHighSchool_bin[JaponW2$ses_education_num < 4] <- 1
JaponW2$ses_educHighSchool_bin[JaponW2$ses_education_num >= 4] <- 0
JaponW2$ses_educCollege_bin <- NA
JaponW2$ses_educCollege_bin[JaponW2$ses_education_num == 4 | JaponW2$ses_education_num == 6] <- 1
JaponW2$ses_educCollege_bin[JaponW2$ses_education_num == 5 | JaponW2$ses_education_num < 4] <- 0
JaponW2$ses_educUni_bin <- NA
JaponW2$ses_educUni_bin[JaponW2$ses_education_num == 5] <- 1
JaponW2$ses_educUni_bin[JaponW2$ses_education_num <= 4 | JaponW2$ses_education_num == 6] <- 0
table(JaponW2$ses_educHighSchool_bin, useNA = "always")
table(JaponW2$ses_educCollege_bin, useNA = "always")
table(JaponW2$ses_educUni_bin, useNA = "always")

### 1.5.5 Personality ####
table(JaponW2$q7_1, useNA = "always")
table(JaponW2$q7_2, useNA = "always")
table(JaponW2$q7_3, useNA = "always")
table(JaponW2$q7_4, useNA = "always")
table(JaponW2$q8, useNA = "always")
JaponW2$ses_extrovertedEnthusiastic_num <- NA
JaponW2$ses_extrovertedEnthusiastic_num[JaponW2$q7_1 == 1] <- 1
JaponW2$ses_extrovertedEnthusiastic_num[JaponW2$q7_1 == 2] <- (2 / 3)
JaponW2$ses_extrovertedEnthusiastic_num[JaponW2$q7_1 == 3] <- (1 / 3)
JaponW2$ses_extrovertedEnthusiastic_num[JaponW2$q7_1 == 4] <- 0
JaponW2$ses_criticalConfrontational_num <- NA
JaponW2$ses_criticalConfrontational_num[JaponW2$q7_2 == 1] <- 1
JaponW2$ses_criticalConfrontational_num[JaponW2$q7_2 == 2] <- (2 / 3)
JaponW2$ses_criticalConfrontational_num[JaponW2$q7_2 == 3] <- (1 / 3)
JaponW2$ses_criticalConfrontational_num[JaponW2$q7_2 == 4] <- 0
JaponW2$ses_reliableDisciplined_num <- NA
JaponW2$ses_reliableDisciplined_num[JaponW2$q7_3 == 1] <- 1
JaponW2$ses_reliableDisciplined_num[JaponW2$q7_3 == 2] <- (2 / 3)
JaponW2$ses_reliableDisciplined_num[JaponW2$q7_3 == 3] <- (1 / 3)
JaponW2$ses_reliableDisciplined_num[JaponW2$q7_3 == 4] <- 0
JaponW2$ses_anxiousIrritable_num <- NA
JaponW2$ses_anxiousIrritable_num[JaponW2$q7_4 == 1] <- 1
JaponW2$ses_anxiousIrritable_num[JaponW2$q7_4 == 2] <- (2 / 3)
JaponW2$ses_anxiousIrritable_num[JaponW2$q7_4 == 3] <- (1 / 3)
JaponW2$ses_anxiousIrritable_num[JaponW2$q7_4 == 4] <- 0
JaponW2$ses_takesRisks_num <- NA
JaponW2$ses_takesRisks_num <- JaponW2$q8 - 1 # reverse code so that higher score means more risk taking
table(JaponW2$ses_extrovertedEnthusiastic_num, useNA = "always")
table(JaponW2$ses_criticalConfrontational_num, useNA = "always")
table(JaponW2$ses_reliableDisciplined_num, useNA = "always")
table(JaponW2$ses_anxiousIrritable_num, useNA = "always")
table(JaponW2$ses_takesRisks_num, useNA = "always")

### 1.5.6 Income ####
table(JaponW2$q26, useNA = "always")
JaponW2$ses_income_num <- JaponW2$q26
JaponW2$ses_income_num[JaponW2$q26 == 16] <- NA
JaponW2$ses_incomeLow_bin[JaponW2$ses_income_num <= 5] <- 1
JaponW2$ses_incomeLow_bin[JaponW2$ses_income_num > 5] <- 0
JaponW2$ses_incomeMid_bin[JaponW2$ses_income_num <= 10 & JaponW2$ses_income_num > 5] <- 1
JaponW2$ses_incomeMid_bin[JaponW2$ses_income_num > 10 | JaponW2$ses_income_num <= 5] <- 0
JaponW2$ses_incomeHigh_bin[JaponW2$ses_income_num > 10] <- 1
JaponW2$ses_incomeHigh_bin[JaponW2$ses_income_num <= 10] <- 0
table(JaponW2$ses_income_num, useNA = "always")
table(JaponW2$ses_incomeLow_bin, useNA = "always")
table(JaponW2$ses_incomeMid_bin, useNA = "always")
table(JaponW2$ses_incomeHigh_bin, useNA = "always")

### 1.5.7 Religion ####
table(JaponW2$q21, useNA = "always")
JaponW2$ses_nonReligious_bin <- 0
JaponW2$ses_nonReligious_bin[JaponW2$q21 == 10] <- 1
JaponW2$ses_nonReligious_bin[JaponW2$q21 == 11 | is.na(JaponW2$q21)] <- NA
table(JaponW2$ses_nonReligious_bin, useNA = "always")

# 2. Correlations between variables of interest ####
## 2.1 AI all ####
JapanBaroAIVariables <- c(
  "q9_1", "q9_2", "q9_3", "q9_4", "q9_5", "q9_6", "q9_7",
  "q10_1", "q10_2", "q10_3", "q10_4", "q10_5", "q10_6", "q10_7",
  "q13_1", "q13_2", "q13_3", "q13_4"
)
modelsummary::datasummary_correlation(JaponW2[JapanBaroAIVariables])
# parmi les questions sur l'IA, les question 9_ corrèlent très bien ensemble et avec q10_1, q10_5, q10_7 eet q13_1
# q13_2 ne matche bien avec rien
JaponW2$dv_technophilia_num <- JaponW2$q9_1 + JaponW2$q9_2 + JaponW2$q9_3 + JaponW2$q9_4 + JaponW2$q9_5 + JaponW2$q9_6 + JaponW2$q9_7 +
  JaponW2$q10_1 + JaponW2$q10_2 + JaponW2$q10_3 + JaponW2$q10_4 + JaponW2$q10_5 + JaponW2$q10_6 + JaponW2$q10_7 +
  JaponW2$q13_1 + JaponW2$q13_2 + JaponW2$q13_3 + JaponW2$q13_4
table(JaponW2$dv_technophilia_num, useNA = "always")

## 2.2 Left-right all ####
JapanBaroLRVariables <- c("q11_5", "q15_1", "q15_2", "q15_3", "q15_4", "leftRight_selfReport_num")
modelsummary::datasummary_correlation(JaponW2[JapanBaroLRVariables])
# pas génial, à part q15_2 et q15_4 ensemble qui loadent très bien (.46) et un peu q15_1 et q15_2 (.22)
# leftRight_selfReport_num loade négativement avec q11_5 et q15_3, ce qui est normal... Mais ces deux-là loadent positivement avec le reste, ce qui est bizarre.

## 2.3 AI dread and controllability
## 2.3.1 AI dread ####
JaponW2$dv_dread_num <- 1 - ((((JaponW2$q9_3 + JaponW2$q9_5 + JaponW2$q9_6 + JaponW2$q9_7) / 4) - 1) / 3)
# reverse-coded so that "Strongly agree" should now be the highest option
AIDread <- c("q9_3", "q9_5", "q9_6", "q9_7")
modelsummary::datasummary_correlation(JaponW2[AIDread])

## 2.3.2 AI controllability ####
JaponW2$dv_controllability_num <- 1 - ((((JaponW2$q10_2 + JaponW2$q10_3 + JaponW2$q10_4 + JaponW2$q10_6) / 4) - 1) / 3)
# reverse-coded so that "Strongly agree" should now be the highest option
AIControllability <- c("q10_2", "q10_3", "q10_4", "q10_6")
modelsummary::datasummary_correlation(JaponW2[AIControllability])

## 2.4 Left-right economics ####
JaponW2$leftRight_3items_num <- 1 - ((((JaponW2$q15_1 + JaponW2$q15_2 + JaponW2$q15_4) / 3) - 1) / 3)
cor.test(JaponW2$leftRight_selfReport_num, JaponW2$leftRight_3items_num) # p<0.001, cor = 0.13
JapanBaroLREconVariables <- c("q15_1", "q15_2", "q15_4")

## 2.5 AI & left-right ####
modelsummary::datasummary_correlation(JaponW2[c(JapanBaroLREconVariables, JapanBaroAIVariables)])
cor.test(JaponW2$leftRight_selfReport_num, JaponW2$dv_dread_num) # N.S., cor = -0.05
cor.test(JaponW2$leftRight_3items_num, JaponW2$dv_dread_num) # p<0.001, cor = 0.21. The more economically left-wing, the higher the dread of AI
cor.test(JaponW2$leftRight_selfReport_num, JaponW2$dv_controllability_num) # N.S., cor = 0.004
cor.test(JaponW2$leftRight_3items_num, JaponW2$dv_controllability_num) # p<0.001, cor = 0.13. The more economically left-wing, the higher the controllability of AI

# 3. Factor analysis ####
## 3.1 Dread factor analysis ####
AIDreadFA <- na.omit(JaponW2[AIDread]) # Retirer les valeurs manquantes
factorAnal <- factanal(AIDreadFA, factors = 1) # Analyse avec 1 facteur
factorLoadingsDread <- as.numeric(factorAnal$loadings[, 1]) # Charges factorielles
(factor1stEigen <- round(eigen(cor(AIDreadFA))$values[1], 2))
# Première valeur propre supérieure à 1, c'est bon!
# Vérification de la cohérence interne avec Cronbach's alpha
(cronbachAlpha <- round(psych::alpha(AIDreadFA)$total$raw_alpha, 2))
# std.alpha de 0.86, pile dans la zone idéale de 0,7 à 0,9!
ggplot(data.frame(AIDread, factorLoadingsDread), aes(x = AIDread, y = factorLoadingsDread)) +
  geom_bar(stat = "identity", fill = "black", color = "black", width = 0.6) +
  coord_flip() +
  geom_hline(yintercept = 0.3, color = "red", linetype = "dashed") +
  geom_text(aes(label = round(factorLoadingsDread, 2)), hjust = -0.2) +
  scale_y_continuous(limits = c(0, 1)) +
  scale_x_discrete(labels = c("Fear AI attack", "Fear AI replace humans", "AI threat to civilization", "Feat AI catastrophe")) +
  annotate("text", label = paste("Cronbach's alpha =", cronbachAlpha), x = 1.55, y = 0.725, size = 4) + # x = 4.5, y = 0.08, size = 6) +
  annotate("text", label = paste("First eigenvalue =", factor1stEigen), x = 0.55, y = 0.725, size = 4) + # x = 4.3, y = 0.08, size = 6) +
  labs(
    title = "Variables associated\nwith each factor",
    x = "Variable",
    y = "Factor loadings"
  ) +
  clessnverse::theme_clean_light()
ggsave("_SharedFolder_article_technophobie/graph/DreadFactorLoad.pdf", width = 5.5, height = 4.25)

## 3.2 Controllability factor analysis ####
AIControllabilityFA <- na.omit(JaponW2[AIControllability]) # Retirer les valeurs manquantes
factorAnal <- factanal(AIControllabilityFA, factors = 1) # Analyse avec 1 facteur
factorLoadingsControllability <- as.numeric(factorAnal$loadings[, 1]) # Charges factorielles
(factor1stEigen <- round(eigen(cor(AIControllabilityFA))$values[1], 2))
# Première valeur propre supérieure à 1, c'est bon!
# Vérification de la cohérence interne avec Cronbach's alpha
(cronbachAlpha <- round(psych::alpha(AIControllabilityFA)$total$raw_alpha, 2))
# std.alpha de 0.7
ggplot(data.frame(AIControllability, factorLoadingsControllability), aes(x = AIControllability, y = factorLoadingsControllability)) +
  geom_bar(stat = "identity", fill = "black", color = "black", width = 0.6) +
  coord_flip() +
  geom_hline(yintercept = 0.3, color = "red", linetype = "dashed") +
  geom_text(aes(label = round(factorLoadingsControllability, 2)), hjust = -0.2) +
  scale_y_continuous(limits = c(0, 1)) +
  scale_x_discrete(labels = c("Experts can control AI", "Predict AI job impact", "Understand AI societal change", "Limit AI risks")) +
  annotate("text", label = paste("Cronbach's alpha =", cronbachAlpha), x = 1.55, y = 0.725, size = 4) + # x = 4.5, y = 0.08, size = 6) +
  annotate("text", label = paste("First eigenvalue =", factor1stEigen), x = 0.55, y = 0.725, size = 4) + # x = 4.3, y = 0.08, size = 6) +
  labs(
    title = "Variables associated\nwith each factor",
    x = "Variable",
    y = "Factor loadings"
  ) +
  clessnverse::theme_clean_light()
ggsave("_SharedFolder_article_technophobie/graph/ControllabilityFactorLoad.pdf", width = 5.5, height = 4.25)

## 3.3 Left-right factor analysis ####
DataForFactAnal <- na.omit(JaponW2[JapanBaroLREconVariables]) # Retirer les valeurs manquantes
factorAnal <- factanal(DataForFactAnal, factors = 1) # Analyse avec 1 facteur
factorLoadings <- as.numeric(factorAnal$loadings[, 1]) # Charges factorielles
(factor1stEigen <- round(eigen(cor(DataForFactAnal))$values[1], 2))
# Première valeur propre supérieure à 1, c'est bon!
# Vérification de la cohérence interne avec Cronbach's alpha
(cronbachAlpha <- round(psych::alpha(DataForFactAnal)$total$raw_alpha, 2))
# std.alpha de 0.53, en-dessous de la zone idéale de 0,7 à 0,9
ggplot(data.frame(JapanBaroLREconVariables, factorLoadings), aes(x = JapanBaroLREconVariables, y = factorLoadings)) +
  geom_bar(stat = "identity", fill = "black", color = "black", width = 0.6) +
  coord_flip() +
  geom_hline(yintercept = 0.3, color = "red", linetype = "dashed") +
  geom_text(aes(label = round(factorLoadings, 2)), hjust = -0.2) +
  scale_y_continuous(limits = c(0, 1)) +
  scale_x_discrete(labels = c("Taxing companies", "Unemployment benefits", "Basic income")) +
  annotate("text", label = paste("Cronbach's alpha =", cronbachAlpha), x = 0.85, y = 0.725, size = 4) +
  annotate("text", label = paste("First eigenvalue =", factor1stEigen), x = 0.65, y = 0.725, size = 4) +
  labs(
    title = NULL,
    x = NULL,
    y = NULL
  ) +
  clessnverse::theme_clean_light()
ggsave("_SharedFolder_article_technophobie/graph/W2LeftRightFactorLoad.pdf", width = 5.5, height = 4.25)

## 3.4 Left-right + AI factor analysis ####

create_factor_loading_table <- function(efa_result, cutoff = 0.30) {
  # Extraire les loadings
  loadings_matrix <- efa_result$loadings[]

  # Convertir en data frame
  loadings_df <- as.data.frame(loadings_matrix)

  # Arrondir à 2 décimales
  loadings_df <- round(loadings_df, 2)

  # Remplacer les valeurs en dessous du seuil par ""
  for (i in 1:ncol(loadings_df)) {
    loadings_df[abs(loadings_df[, i]) < cutoff, i] <- ""
  }

  # Traitement spécial pour q15_1 : forcer les cellules vides dans Factor 1 et 2
  if ("q15_1" %in% rownames(loadings_df)) {
    # Sauvegarder la valeur originale de Factor 3
    q15_1_factor3 <- loadings_matrix["q15_1", 3]

    # Vider Factor 1 et Factor 2 pour q15_1
    loadings_df["q15_1", 1] <- "" # Factor 1
    loadings_df["q15_1", 2] <- "" # Factor 2

    # Garder la valeur de Factor 3 (même si < 0.30)
    loadings_df["q15_1", 3] <- round(q15_1_factor3, 2)
  }

  # Créer les libellés des questions
  question_labels <- c(
    "q10_2" = "Experts can control AI",
    "q10_3" = "Predict AI job impact",
    "q10_4" = "Understand AI societal change",
    "q10_6" = "Limit AI risks",
    "q9_3" = "Fear AI attack",
    "q9_5" = "Fear AI replace humans",
    "q9_6" = "AI threat to civilization",
    "q9_7" = "Fear AI catastrophe",
    "q15_1" = "Taxing companies",
    "q15_2" = "Unemployment benefits",
    "q15_4" = "Basic income"
  )

  # Ajouter les descriptions
  loadings_df$`Item Description` <- question_labels[rownames(loadings_df)]

  # Réorganiser les colonnes (mettre la description en premier)
  loadings_df <- loadings_df[, c(ncol(loadings_df), 1:(ncol(loadings_df) - 1))]

  # Renommer les colonnes des facteurs
  colnames(loadings_df)[2:ncol(loadings_df)] <- paste("Factor", 1:(ncol(loadings_df) - 1))

  # Ordonner les lignes selon l'ordre souhaité
  desired_order <- c("q10_2", "q10_3", "q10_4", "q10_6", "q9_3", "q9_5", "q9_6", "q9_7", "q15_1", "q15_2", "q15_4")
  existing_order <- desired_order[desired_order %in% rownames(loadings_df)]
  loadings_df <- loadings_df[existing_order, ]

  return(loadings_df)
}

DataForFactAnal <- na.omit(JaponW2[c(AIControllability, AIDread, JapanBaroLREconVariables)]) # Retirer les valeurs manquantes


# Extraire plusieurs solutions avec différents nombres de facteurs

efa_1 <- fa(DataForFactAnal, nfactors = 1, rotate = "oblimin", fm = "pa")
efa_2 <- fa(DataForFactAnal, nfactors = 2, rotate = "oblimin", fm = "pa")
efa_3 <- fa(DataForFactAnal, nfactors = 3, rotate = "oblimin", fm = "pa")
efa_4 <- fa(DataForFactAnal, nfactors = 4, rotate = "oblimin", fm = "pa")
# efa_5 <- fa(DataForFactAnal, nfactors = 5, rotate = "oblimin", fm = "pa")

(factor1stEigen <- round(eigen(cor(DataForFactAnal))$values[1], 2))


# Première valeur propre supérieure à 1, c'est bon!
# Vérification de la cohérence interne avec Cronbach's alpha
(cronbachAlpha <- round(psych::alpha(DataForFactAnal)$total$raw_alpha, 2))

# std.alpha de 0.68, un peu en dessous de la zone idéale de 0,7 à 0,9

# Comparer les indices d'ajustement

rbind(
  "1 facteur" = c(RMSEA = efa_1$RMSEA[1], TLI = efa_1$TLI, variance = sum(efa_1$Vaccounted[2, 1:1])),
  "2 facteurs" = c(RMSEA = efa_2$RMSEA[1], TLI = efa_2$TLI, variance = sum(efa_2$Vaccounted[2, 1:2])),
  "3 facteurs" = c(RMSEA = efa_3$RMSEA[1], TLI = efa_3$TLI, variance = sum(efa_3$Vaccounted[2, 1:3])),
  "4 facteurs" = c(RMSEA = efa_4$RMSEA[1], TLI = efa_4$TLI, variance = sum(efa_4$Vaccounted[2, 1:4]))
)

# Créer le tableau
factor_table <- create_factor_loading_table(efa_3, cutoff = 0.30)

# Afficher avec kable
q15_1_row <- which(rownames(factor_table) == "q15_1")

# Afficher avec kable
tableau_latex <- knitr::kable(factor_table,
  format = "latex", # Spécifier le format LaTeX
  digits = 2,
  caption = "Table A2. Exploratory Factor Analysis (Japan)",
  align = c("l", rep("c", ncol(factor_table) - 1)),
  escape = FALSE,
  booktabs = TRUE
) |> # Pour un meilleur style LaTeX
  kableExtra::kable_styling(latex_options = c("striped", "hold_position")) |>
  kableExtra::row_spec(0, bold = TRUE) |>
  kableExtra::row_spec(q15_1_row, bold = TRUE) |>
  kableExtra::column_spec(1, width = "6cm")

knitr::kable(factor_table,
  format = "latex",
  digits = 2,
  booktabs = TRUE,
  escape = FALSE
)



# 4. Scatterplots Left-right + AI ####
# Function to get ribbon range from geom_smooth
get_ribbon_range <- function(data, y_var, x_var, method) {
  temp_data <- data[complete.cases(data[, c(y_var, x_var)]), ]

  if (method == "lm") {
    model <- lm(as.formula(paste(y_var, "~", x_var)), data = temp_data)
    pred <- predict(model, interval = "confidence")
    return(range(pred[, c("lwr", "upr")]))
  } else if (method == "loess") {
    model <- loess(as.formula(paste(y_var, "~", x_var)), data = temp_data)
    pred <- predict(model, se = TRUE)
    ci_lower <- pred$fit - 1.96 * pred$se.fit
    ci_upper <- pred$fit + 1.96 * pred$se.fit
    return(range(c(ci_lower, ci_upper), na.rm = TRUE))
  }
}

# Get ranges from all ribbons
ranges <- c(
  get_ribbon_range(CanW2, "dv_dread_num", "leftRight_selfReport_num", "lm"),
  get_ribbon_range(CanW2, "dv_dread_num", "leftRight_selfReport_num", "loess"),
  get_ribbon_range(CanW2, "dv_dread_num", "leftRight_3items_num", "lm"),
  get_ribbon_range(CanW2, "dv_dread_num", "leftRight_3items_num", "loess"),
  get_ribbon_range(CanW2, "dv_controllability_num", "leftRight_selfReport_num", "lm"),
  get_ribbon_range(CanW2, "dv_controllability_num", "leftRight_selfReport_num", "loess"),
  get_ribbon_range(CanW2, "dv_controllability_num", "leftRight_3items_num", "lm"),
  get_ribbon_range(CanW2, "dv_controllability_num", "leftRight_3items_num", "loess")
)

y_range <- range(ranges, na.rm = TRUE)

p1 <- ggplot(JaponW2, aes(y = dv_dread_num, x = leftRight_selfReport_num)) +
  # geom_point(position = "jitter", size = 0.07) +
  geom_smooth(method = "lm", color = "#6e6e6d") +
  geom_smooth(method = "loess", color = "black") +
  labs(title = "", y = "AI Dread", x = "Ideology: Right-Leaning\n(Self-Report)") +
  clessnverse::theme_clean_light() +
  theme(text = element_text(family = "serif"))

p2 <- ggplot(JaponW2, aes(y = dv_dread_num, x = leftRight_3items_num)) +
  # geom_point(position = "jitter", size = 0.07) +
  geom_smooth(method = "lm", color = "#6e6e6d") +
  geom_smooth(method = "loess", color = "black") +
  labs(title = "", y = "AI Dread", x = "Ideology: Fiscally Conservative\n(3-Item Scale)") +
  clessnverse::theme_clean_light() +
  theme(text = element_text(family = "serif"))

p3 <- ggplot(JaponW2, aes(y = dv_controllability_num, x = leftRight_selfReport_num)) +
  # geom_point(position = "jitter", size = 0.07) +
  geom_smooth(method = "lm", color = "#6e6e6d") +
  geom_smooth(method = "loess", color = "black") +
  labs(title = "", y = "AI Controllability", x = "Ideology: Right-Leaning\n(Self-Report)") +
  clessnverse::theme_clean_light() +
  theme(text = element_text(family = "serif"))

p4 <- ggplot(JaponW2, aes(y = dv_controllability_num, x = leftRight_3items_num)) +
  # geom_point(position = "jitter", size = 0.07) +
  geom_smooth(method = "lm", color = "#6e6e6d") +
  geom_smooth(method = "loess", color = "black") +
  labs(title = "", y = "AI Controllability", x = "Ideology: Fiscally Conservative\n(3-Item Scale)") +
  clessnverse::theme_clean_light() +
  theme(text = element_text(family = "serif"))

cor_plot <- p1 + p2 + p3 + p4 +
  plot_layout(nrow = 2)
cor_plot
ggsave("_SharedFolder_article_technophobie/graph/patchwork_plot.pdf", plot = cor_plot, width = 6, height = 6, dpi = 600)
ggsave("_SharedFolder_article_technophobie/graph/patchwork_plot.jpeg", plot = cor_plot, width = 6, height = 6, dpi = 600)

# 5. Regressions ####
## 5.1 Dread ####
ModDreadSelfRep <- lm(data = JaponW2, dv_dread_num ~ leftRight_selfReport_num)
summary(ModDreadSelfRep) # N.S., -0.05
ModDreadCreated <- lm(data = JaponW2, dv_dread_num ~ leftRight_3items_num)
summary(ModDreadCreated) # p<0.001, 0.23
ModDreadSelfRepSES <- lm(data = JaponW2, dv_dread_num ~ leftRight_selfReport_num + ses_female_bin + ses_age3555_bin + ses_age55m_bin + ses_incomeMid_bin + ses_incomeHigh_bin +
  ses_educCollege_bin + ses_educUni_bin)
summary(ModDreadSelfRepSES)
ModDreadCreatedSES <- lm(data = JaponW2, dv_dread_num ~ leftRight_3items_num + ses_female_bin + ses_age3555_bin + ses_age55m_bin + ses_incomeMid_bin + ses_incomeHigh_bin +
  ses_educCollege_bin + ses_educUni_bin)
summary(ModDreadCreatedSES)
ModDreadSelfRepCtrl <- lm(data = JaponW2, dv_dread_num ~ leftRight_selfReport_num + ses_female_bin + ses_age3555_bin + ses_age55m_bin + ses_incomeMid_bin + ses_incomeHigh_bin +
  ses_educCollege_bin + ses_educUni_bin + ses_takesRisks_num + ses_extrovertedEnthusiastic_num + ses_criticalConfrontational_num +
  ses_reliableDisciplined_num + ses_anxiousIrritable_num)
summary(ModDreadSelfRepCtrl) # N.S., -0.04
ModDreadCreatedCtrl <- lm(data = JaponW2, dv_dread_num ~ leftRight_3items_num + ses_female_bin + ses_age3555_bin + ses_age55m_bin + ses_incomeMid_bin + ses_incomeHigh_bin +
  ses_educCollege_bin + ses_educUni_bin + ses_takesRisks_num + ses_extrovertedEnthusiastic_num + ses_criticalConfrontational_num +
  ses_reliableDisciplined_num + ses_anxiousIrritable_num)
summary(ModDreadCreatedCtrl) # p<0.001, 0.17
modelsummary::modelsummary(
  list(
    ModDreadSelfRep, ModDreadCreated,
    ModDreadSelfRepSES, ModDreadCreatedSES, ModDreadSelfRepCtrl, ModDreadCreatedCtrl
  ),
  statistic = "std.error", stars = TRUE,
  coef_rename = c(
    "leftRight_selfReport_num" = "Ideology: Right-Leaning (Self-Report)",
    "leftRight_3items_num" = "Ideology: Fiscally Conservative (3-Item Scale)",
    "ses_female_bin" = "Gender: Woman",
    "ses_age3555_bin" = "Age: 35-54",
    "ses_age55m_bin" = "Age: 55+",
    "ses_incomeMid_bin" = "Income: Middle",
    "ses_incomeHigh_bin" = "Income: High",
    "ses_educCollege_bin" = "Education: College",
    "ses_educUni_bin" = "Education: University",
    "ses_extrovertedEnthusiastic_num" = "Personality: Extroverted, Enthusiastic",
    "ses_takesRisks_num" = "Takes Risks",
    "ses_criticalConfrontational_num" = "Personality: Critical, Confrontational",
    "ses_reliableDisciplined_num" = "Personality: Reliable, Disciplined",
    "ses_anxiousIrritable_num" = "Personality: Anxious, Irritable"
  ),
  title = "Variables influencing AI dread in Japan \\label{table:DreadVariablesJap}", escape = FALSE, gof_omit = "Log.Lik.|RMSE",
  notes = c("OLS regression. \\textit{Reference categories}: Gender: Man, Age: 18-34, Income: Low, Education: High School or Less, Does Not Take Risks, Personality: Not [characteristic]"),
  output = "latex"
)
summary(lm(data = JaponW2, dv_dread_num ~ leftRight_selfReport_num_no5)) # still non-significant
summary(lm(data = JaponW2, dv_dread_num ~ leftRight_selfReport_num_no5 + ses_female_bin + ses_age3555_bin + ses_age55m_bin + ses_incomeMid_bin + ses_incomeHigh_bin +
  ses_educCollege_bin + ses_educUni_bin + ses_takesRisks_num + ses_extrovertedEnthusiastic_num + ses_criticalConfrontational_num +
  ses_reliableDisciplined_num + ses_anxiousIrritable_num)) # still non-significant

# 5.2 Controllability ####
ModControllabilitySelfRep <- lm(data = JaponW2, dv_controllability_num ~ leftRight_selfReport_num)
summary(ModControllabilitySelfRep) # N.S., 0.003
ModControllabilityCreated <- lm(data = JaponW2, dv_controllability_num ~ leftRight_3items_num)
summary(ModControllabilityCreated) # p<0.001, 0.11
ModControllabilitySelfRepSES <- lm(data = JaponW2, dv_controllability_num ~ leftRight_selfReport_num + ses_female_bin + ses_age3555_bin + ses_age55m_bin + ses_incomeMid_bin + ses_incomeHigh_bin +
  ses_educCollege_bin + ses_educUni_bin)
summary(ModControllabilitySelfRepSES)
ModControllabilityCreatedSES <- lm(data = JaponW2, dv_controllability_num ~ leftRight_3items_num + ses_female_bin + ses_age3555_bin + ses_age55m_bin + ses_incomeMid_bin + ses_incomeHigh_bin +
  ses_educCollege_bin + ses_educUni_bin)
summary(ModControllabilityCreatedSES)
ModControllabilitySelfRepCtrl <- lm(data = JaponW2, dv_controllability_num ~ leftRight_selfReport_num + ses_female_bin + ses_age3555_bin + ses_age55m_bin + ses_incomeMid_bin + ses_incomeHigh_bin +
  ses_educCollege_bin + ses_educUni_bin + ses_takesRisks_num + ses_extrovertedEnthusiastic_num + ses_criticalConfrontational_num +
  ses_reliableDisciplined_num + ses_anxiousIrritable_num)
summary(ModControllabilitySelfRepCtrl) # N.S., 0.008
ModControllabilityCreatedCtrl <- lm(data = JaponW2, dv_controllability_num ~ leftRight_3items_num + ses_female_bin + ses_age3555_bin + ses_age55m_bin + ses_incomeMid_bin + ses_incomeHigh_bin +
  ses_educCollege_bin + ses_educUni_bin + ses_takesRisks_num + ses_extrovertedEnthusiastic_num + ses_criticalConfrontational_num +
  ses_reliableDisciplined_num + ses_anxiousIrritable_num)
summary(ModControllabilityCreatedCtrl) # p<0.05, 0.08
modelsummary::modelsummary(
  list(
    ModControllabilitySelfRep, ModControllabilityCreated,
    ModControllabilitySelfRepSES, ModControllabilityCreatedSES,
    ModControllabilitySelfRepCtrl, ModControllabilityCreatedCtrl
  ),
  statistic = "std.error", stars = TRUE,
  coef_rename = c(
    "leftRight_selfReport_num" = "Ideology: Right-Leaning (Self-Report)",
    "leftRight_3items_num" = "Ideology: Fiscally Conservative (3-Item Scale)",
    "ses_female_bin" = "Gender: Woman",
    "ses_age3555_bin" = "Age: 35-54",
    "ses_age55m_bin" = "Age: 55+",
    "ses_incomeMid_bin" = "Income: Middle",
    "ses_incomeHigh_bin" = "Income: High",
    "ses_educCollege_bin" = "Education: College",
    "ses_educUni_bin" = "Education: University",
    "ses_extrovertedEnthusiastic_num" = "Personality: Extroverted, Enthusiastic",
    "ses_takesRisks_num" = "Takes Risks",
    "ses_criticalConfrontational_num" = "Personality: Critical, Confrontational",
    "ses_reliableDisciplined_num" = "Personality: Reliable, Disciplined",
    "ses_anxiousIrritable_num" = "Personality: Anxious, Irritable"
  ),
  title = "Variables influencing AI controllability in Japan \\label{table:ControllabilityVariablesJap}", escape = FALSE, gof_omit = "Log.Lik.|RMSE",
  notes = c("OLS regression. \\textit{Reference categories}: Gender: Man, Age: 18-34, Income: Low, Education: High School or Less, Does Not Take Risks, Personality: Not [characteristic]"),
  output = "latex"
)
summary(lm(data = JaponW2, dv_controllability_num ~ leftRight_selfReport_num_no5)) # still non-significant
summary(lm(data = JaponW2, dv_controllability_num ~ leftRight_selfReport_num_no5 + ses_female_bin + ses_age3555_bin + ses_age55m_bin + ses_incomeMid_bin + ses_incomeHigh_bin +
  ses_educCollege_bin + ses_educUni_bin + ses_takesRisks_num + ses_extrovertedEnthusiastic_num + ses_criticalConfrontational_num +
  ses_reliableDisciplined_num + ses_anxiousIrritable_num)) # still non-significant

# ===================== 6. Predicted probabilities ========================
## 6.1 Dread ####
ModDreadSelfRepCtrl
ModDreadCreatedCtrl

# load CM Roman font
extrafont::font_import()
extrafont::font_install("fontcm")
extrafont::loadfonts()

dread_selfReport_pred <- marginaleffects::plot_predictions(ModDreadSelfRepCtrl, condition = "leftRight_selfReport_num") +
  labs(
    x = "Left-Right Self-Report",
    y = "Dread Scale"
  ) +
  clessnverse::theme_clean_light() +
  theme(text = element_text(family = "CM Roman"))
dread_selfReport_pred
ggsave("_SharedFolder_article_technophobie/graph/dread_selfReport_pred.pdf", plot = dread_selfReport_pred, width = 5.5, height = 4.25)

dread_Created_pred <- marginaleffects::plot_predictions(ModDreadCreatedCtrl, condition = "leftRight_3items_num") +
  labs(
    x = "Fiscal Conservatism",
    y = "Dread Scale"
  ) +
  clessnverse::theme_clean_light() +
  theme(text = element_text(family = "CM Roman"))
dread_Created_pred
ggsave("_SharedFolder_article_technophobie/graph/dread_Created_pred.pdf", plot = dread_Created_pred, width = 5.5, height = 4.25)

## 6.2 Controllability ####
ModControllabilitySelfRepCtrl
ModControllabilityCreatedCtrl

Controllability_selfReport_pred <- plot_predictions(ModControllabilitySelfRepCtrl, condition = "leftRight_selfReport_num") +
  labs(
    x = "Left-Right Self-Report",
    y = "Controllability Scale"
  ) +
  clessnverse::theme_clean_light() +
  theme(text = element_text(family = "CM Roman"))
Controllability_selfReport_pred
ggsave("_SharedFolder_article_technophobie/graph/Controllability_selfReport_pred.pdf", plot = Controllability_selfReport_pred, width = 5.5, height = 4.25)

Controllability_Created_pred <- plot_predictions(ModControllabilityCreatedCtrl, condition = "leftRight_3items_num") +
  labs(
    x = "Fiscal Conservatism",
    y = "Controllability Scale"
  ) +
  clessnverse::theme_clean_light() +
  theme(text = element_text(family = "CM Roman"))
Controllability_Created_pred
ggsave("_SharedFolder_article_technophobie/graph/Controllability_Created_pred.pdf", plot = Controllability_Created_pred, width = 5.5, height = 4.25)
predicted_plot <- dread_selfReport_pred + dread_Created_pred + Controllability_selfReport_pred + Controllability_Created_pred +
  plot_layout(nrow = 2)
predicted_plot
ggsave("_SharedFolder_article_technophobie/graph/predicted_plot.pdf", plot = predicted_plot, width = 6, height = 6)

# 7. Descriptive stats ####
# Create the plots first
DreadDescriptive <- ggplot(JaponW2, aes(x = dv_dread_num)) +
  geom_histogram(bins = 11, fill = "black", color = "white") +
  labs(x = "AI Dread Score", y = "Frequency") +
  clessnverse::theme_clean_light() +
  theme(text = element_text(family = "serif"))
# do same graph but with y-axis = percentage
# ggplot(JaponW2, aes(x = dv_dread_num)) +
#  geom_histogram(aes(y = after_stat(count) / sum(after_stat(count))), bins = 11, fill = "black", color = "white") +
#  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
#  labs(x = "AI Dread Score", y = "Percentage") +
#  clessnverse::theme_clean_light() +
#  theme(text = element_text(family = "CM Roman")) +
#  coord_cartesian(ylim = c(0, max_y))

ControllabilityDescriptive <- ggplot(JaponW2, aes(x = dv_controllability_num)) +
  geom_histogram(bins = 11, fill = "black", color = "white") +
  labs(x = "AI Controllability Score", y = "Frequency") +
  clessnverse::theme_clean_light() +
  theme(text = element_text(family = "serif"))
FiscalConservatismDescriptive <- ggplot(JaponW2, aes(x = leftRight_3items_num)) +
  geom_histogram(bins = 10, fill = "black", color = "white") +
  labs(x = "Fiscal Conservatism Score", y = "Frequency") +
  clessnverse::theme_clean_light() +
  theme(text = element_text(family = "serif"))
LeftRightDescriptive <- ggplot(JaponW2, aes(x = leftRight_selfReport_num)) +
  geom_histogram(bins = 10, fill = "black", color = "white") +
  labs(x = "Left-Right Self-Report Score", y = "Frequency") +
  clessnverse::theme_clean_light() +
  theme(text = element_text(family = "serif"))

# Find max y value across all vars
max_y <- max(
  ggplot_build(DreadDescriptive)$data[[1]]$count,
  ggplot_build(ControllabilityDescriptive)$data[[1]]$count,
  ggplot_build(FiscalConservatismDescriptive)$data[[1]]$count,
  ggplot_build(LeftRightDescriptive)$data[[1]]$count
)

# Update plots with consistent y-axis
DreadDescriptive <- DreadDescriptive + coord_cartesian(ylim = c(0, max_y))
ControllabilityDescriptive <- ControllabilityDescriptive + coord_cartesian(ylim = c(0, max_y))
FiscalConservatismDescriptive <- FiscalConservatismDescriptive + coord_cartesian(ylim = c(0, max_y))
LeftRightDescriptive <- LeftRightDescriptive + coord_cartesian(ylim = c(0, max_y))
# combine all four graphs in a 2x2 grid
descriptive_stats_plot <- DreadDescriptive + ControllabilityDescriptive + FiscalConservatismDescriptive + LeftRightDescriptive +
  plot_layout(nrow = 2)
descriptive_stats_plot
ggsave("_SharedFolder_article_technophobie/graph/descriptive_stats_plot.pdf", plot = descriptive_stats_plot, width = 6, height = 6)

prop.table(table(JaponW2$leftRight_selfReport_num, useNA = "always"))
mean(JaponW2$dv_dread_num, na.rm = TRUE) # 4.04
mean(JaponW2$dv_controllability_num, na.rm = TRUE) # 4.04
mean(JaponW2$leftRight_selfReport_num, na.rm = TRUE) # 4.04
mean(JaponW2$leftRight_3items_num, na.rm = TRUE) # 4.04
